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■ Abstract 

Observations show that the magnetic field in young supernova remnants (SNRs) is significantly stronger than can be 
expected fi'om the compression of the circumstellar medium (CSM) by a factor of four expected for strong blast waves. 
Additionally, the polarization is mainly radial, which is also contrary to expectation from compression of the CSM 
magnetic field. Cosmic rays (CRs) may help to explain these two observed features. They can increase the compression 
ratio to factors well over those of regular strong shocks by adding a relativistic plasma component to the pressure, and 
|by draining the shock of energy when CRs escape from the region. The higher compression ratio will also allow for the 
contact discontinuity, which is subject to the Rayleigh- Taylor (R-T) instability, to reach much further out to the forward 
shock. This could create a preferred radial polarization of the magnetic field. With an adaptive mesh refinement MHD 
code (AMRVAC), we simulate the evolution of SNRs with three different configurations of the initial CSM magnetic 
field, and look at two different equations of state in order to look at the possible influence of a CR plasma component. 
The spectrum of CRs can be simulated using test particles, of which we also show some preliminary results that agree 
well with available analytical solutions. 
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1. Introduction 

Supernova remnants (SNR) are sites where particles 
are believed to be accelerated up to energies of 10^^ eV. 
.Cosmic rays (CRs) scatter off local magnetic field per- 
turbations and are accelerated by the first order Fermi 
.mechanism, also referred to as diffusive shock accelera- 
tion (DSA). In order to accelerate CRs to sufficiently high 
energies and explain the non-thermal X-ray synchrotron 
emission that is confined to a small region around the 
blast wave, the magnetic field needs to be amplified b y 
a factor of about 100 (e.g. IVolk et all l2005t IVinkl . l2005l ) . 
This is much more than can be expected by mere com- 
pression of the circumstellar magnetic field by a factor 
four in a typical strong blast wave. The spectral proper- 
ties of the non-thermal emission indicate that the cosmic 
ray diffusion must be near the Boh m-hmit a n d that mag- 



netic field turbulence is high (e.g. iVinkl (|2005[ ): IStage et al 
(200(1)). Additionally, information on the field topology 
from radio-synchrotron polarization data indicates that 
young SNRs have a mostly radially oriented magnetic field, 
whereas the field in old r emnants is mainly circumferential 
(jPickel and Milnel . [l976h . 



between CRs and magnetic fields in SNRs, we believe it 
is important to also understand the structure of magnetic 
fields better . Although th e CRs may self-generate mag- 
netic fields ( Skillingl . [l97^ , for now we will treat them as 
test-particles and only look at the evolution of the field 
assuming ideal MHD. 

The contact discontinuity between the supernova ejecta 
and the compressed circumstellar medium, just behind the 
SNR blast wave, is unstable to the Ray l eigh- T aylor (R- 
T) instability. It was suggested by ICullI (|l973l ) that the 
radial fingers that result from this instability stretch the 
field lines radially, thus creating a dominant radial mag- 
netic field component. The stretching of the field along 
the R-T fin gers has been c onfirmed by simulations, for 
instance see Jun et aP (l995) . While this causes the dom- 
inance of the radial field signature in part of the remnant, 
for a typical shock compression ratio of four, it does not 
explain the radial field polarization observed near the for- 
ward shock, since the R-T fingers do not extend that far 



outward ( Chevalier et al 



1992 



Blondin and EllisonI ( 2001 ) have performed hydrody- 



The diffusion of CRs in SNRs depends on the orienta- 
tion of the mean field and the level of magnetic turbulence. 
In order to obtain a better understanding of the interplay 



Email address: K.M.Schure@phys.uu.nl (K.M. Schure) 



namic simulations of SNRs in which they show that, when 
the compression ratio is higher, the R-T fingers almost 
reach out to the blast wave. The presence of a cosmic 
ray gas with sufficient pressure would lead to an effective 
adiabatic index lower than that for a mono-atomic non- 
relativistic gas, which has 7 = 5/3, and thus to a higher 
compression ratio (s = (7 -I- l)/(7 — 1)). For a gas that 
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is dominated by relativistic particles, such as cosmic rays, 
the adiabatic index 7 = 4/3. AdditionaUy, escaping cosmic 
rays drain energy from the shock region, thereby further 
increasing the compression ratio. 

We explore the development of Rayleigh- Taylor fingers 
and magnetic field variations of a SNR inside a stellar wind 
medium with a density p oc r"^ for different initial mag- 
netic topologies. We also investigate the influence of a 
softened equation of state on the properties of the result- 
ing SNR. We compare our results with those from various 
studies that have been performed on this subject in a ho- 
mogeneous and/or unmagnetised medium. The models for 
th e magnetic fiel d are t he same ones we used and described 
by ISchure et al.l ( 2008h . They will be summarised in sec- 
tion [21 where we will also describe the implementation of 
the problem in the code we use: AMRVAC. In line with the 
hydrodynamic simulations by Blondin and Ellison (2001) 
we adopt a softer equation of state, and compare the dif- 
ferences of our MHD models for the values of 7 = 5/3 and 
7 = 1.1, an extreme case that is relevant if cosmic rays are 
the dominant component of the plasma and escape, soft- 
ening the effective equation of state to one with 7 < 4/3. 
Our focus will be on the early evolution of the SNR, in- 
vestigating if we can reproduce the observed radial field 
based on a choice of initial magnetic field topologies and 
the equation of state. 



2. Method 

We use the same model for th e SN ejecta and cir - 
cumstellar wind as that adopted in ISchure et al.l (|2008l ) , 
in which the supernova ejecta evolve in a circumstellar 
medium (GSM) that is shaped by a pre-supernova wind. 
The initial grid is filled with a RSG wind, with a mass 
loss rate of M — 1.54 x 10~^ Mq yr~^, the wind veloc- 
ity V = 4.7 km s^^ (at R w 0.03 pc) and the tempera- 
ture T — 1000 K. The SN ejecta consist of a constant- 
density core, with an envelope for which the density de- 
creases as p oc r^^, which is the typical density profile of 
the ejecta af ter propagation through the star in explosion 
models (c.f. iTruelove and McKeeLll999l ) . The explosion 
mass and energy are respectively set to be: Afoj = 2.5Af0 



The divergence of the magnetic field is controlled by 
adding a source term proportional to V • i? to the induc- 
tion equation, wh ile maintaining conse r vation of momen - 
tum and energy ( Keppens et al. . 2003; Janhunen . 2000l ). 
Additionally, the initial magnetic field strength is chosen 
to be weak, such that it does not infiuence the dynamics. 

The initial grid is filled with a RSG wind, after which 
we introduce a magnetic field into the entire grid, and 
the energetic supernova ejecta into the inner 0.3 pc. We 
implement three different models for the initial magnetic 
field topology, which are described in ISchure et~al.. (2008i) 
and will be summarised below: 

In Models A and D the field is predominantly toroidal, 
with a small radial component, and is modelled after a ro- 
tating dipole field in the equatorial plane with a sin 9 de- 
pcndenc e such tha t the toroidal fiel d vanishe s at the poles 
(Glicvali er and L uo. 1994; Garcia-S egura et a l.. 1999). The 
different components of the field decay with radius as i?0 oc 



^ and Br 



oc r 



such that away from the star, at ra- 



dial scales relevant to the SNR evolution, the toroidal field 
component strongly dominates. 

In Model B and E, we add a 2D turbulent magnetic 
field at the time when we introduce the ejecta. This is 
to mimic the presence of cosmic rays that may induce 
such turbulence. The azimuthal i?^ component is initially 
zero. The field is calculat ed on a 2D cartesian grid t o 
be divergence-free, following lGiacalone and ,Tokipiil|l999l ). 
and transcribed to spherical coordinates in the r, 9 plane. 
The field is set up using the real part of: 



SB{x,y) = ^ Aok~"'^^i{cos(j)ny - sm(j)nx) (5) 

gikn (x COS (f>n~\-y sin (f>n)+il3n 



n=l 



and Eej = 2 x 10^^ erg. The maximum velocity of the 
outer ejecta is chosen to be 15,000 km s~^. 

We use the Adaptive Mesh Re finement version of the 
Versa tile Advection Gode: AMRVAC (' van der Hoist and Keppeni 
20071 ) to solve the MHD equations in the r — 9 plane of a 
spherical grid, with symmetry around the polar axis. The 
equations that are solved are the conservative equations for 
mass, momentum and energy, and the induction equation: 



dtp + \7- iyp) = (1) 

dt (pv) + V • (vpv) - BB) + Vptot = (2) 

ate + V- (ve-BB-v-Hvptot) = (3) 

atB-(- V- (vB-Bv) = -(V-B)v (4) 



X e 

with A^nfe = 256 different modes with wavenumbers k. The 
phase (3 and polarization (j) is randomly chosen between 
and 27r for each wavenumber. The wavenumbers are log- 
arithmically spaced and span the range between a wave- 
length corresponding to two gridcells, and a wavelength 
that covers the entire grid. We approximate a Kolmogorov 
spectrum by adopting a — 5/3. 

In Model G and F, the magnetic field is taken to be 
unidirectional and parallel, representing a situation like 
the dominant ordered magnetic field as observed in e.g. 
%)iral galaxies. We choose the direction of the field to 
be alighned with the symmetry axis. This magnetic field 
topology may in reality be more appropriate to model type 
la SNe instead of the ones that explode in a stellar-wind 
environment. In spherical coordinates, this uniform mag- 
netic field is represented by Br = Bcos9, Bg = —Bsm9. 

Models A-G are set up with an adiabatic index of 7 = 
5/3, the value for a monatomic gas. The grid spans 1.8 x 
10^^ cm radially, and an effective resolution of 2.75 x 
10^^ cm by 0.125° is reached in regions where strong den- 
sity and/or velocity gradients are present. When cosmic 



2 



ray pressure dominates the fluid, the value of the adia- 
batic index drops, approaching the value 7 = 4/3 for a 
relativistically hot gas. Because of escape of cosmic rays 
the equation of state can soften even further. In this pa- 
per we evaluate the extreme case where 7 = 1.1 in models 
D-F. In those simulations, the grid spans 3.0 x 10^^ cm 
radially and the effective resolution reaches 6.25 x 10^^ cm 
by 0.125°. 

3. Results 

The results from simulations of Models A, B and C, are 
plotted in Figures [T] to [31 The density, radial and toroidal 
magnetic field, and additionally the ratio of the radial to 
the total field are plotted for the three assumed initial 
magnetic-field topologies. The results shown here are for 
remnants that have evolved for a period of 634 year. 
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Figure 1: Simulation of SNR into CSM with a mostly toroidal initial 
field (Model A) and an adiabatic index of 7 = 5/3, shown at a 
time of 634 yr after explosion. The upper left quadrant shows the 
logarithm of the density, which shows clearly the occurence of the 
Rayleigh- Taylor instability at the contact discontinuity. The density 
jumps at a radius of 1.3 X 10^^ cm and ~ 7.0 X 10^* cm correspond 
to the locations of the forward and the reverse shock respectively. 
The lower left panel shows the logarithm of the absolute value of the 
radial magnetic field. The lower right part shows the logarithm of 
I I , and the upper right quadrant shows the logarithm of the ratio 
of the radial field component relative to the total field (|i3r|/|B|). 

The supernova ejecta sweeps up the circumstellar mat- 
ter into a dense shell. Four distinct regions can be iden- 
tified: the unshocked circumstellar medium (CSM) ahead 
of the blast wave, the shocked CSM, the shocked ejecta. 




Figure 2: Simulation of the SNR into CSM with a turbulent initial 
field (Model B). Similar to Figure^ the density and magnetic field 
is plotted. Note that in this case the lower right part shows the 
logarithm of \Bg\, and the upper right quadrant shows the ratio of 
the radial field component relative to the total field (|iJ,. |/|B|). 
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Figure 3: Simulation of SNR into CSM with an initial magnetic field 
parallel to the symmetry axis, i.e. vertically aligned (Model C). The 
upper quadrants show again the logarithm of the density and the 
ratio of radial field over the total field. The lower quadrants show 
the absolute values of Br and Bg. 
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and the freely expanding ejecta. The shocked CSM is sep- 
arated from the unshocked CSM by a strong shock (in 
Figures [T][3] at a radius of R[ w 1.3 x 10^^ cm), charac- 
terised by a pressure jump, which increases the density 
by a factor (7 -I- 1)/ (7 — 1) and decreases the velocity in 
the shock frame by the same amount. The shocked CSM is 
separated from the third region, the shocked ejecta, by the 
contact discontinuity, characterised by a jump in the den- 
sity, but at constant pressure. This is the location where 
the R-T instability develops. The reverse shock, at a ra- 
dius of i?r ~ 6.5 X 10^^ cm, marks the boundary with the 
unshocked ejecta. 

The deceleration of the shocked ejecta by the less-dense 
shocked CSM is Rayleigh- Taylor unstable. Since the mag- 
netic field is very weak, it does not influence the dynamics 
of the blast wave, and we do not sec a difference between 
the three models in the development of the R-T instabil- 
ities and the propagation of the blast wave. The com- 
pression ratio at the forward shock is too low to allow 
for the Rayleigh- Taylor fingers to extend out to the for 



ward sh ock, in agreement with finding s bv lChevalier et al 
Chevalier and BlondinI (|199 5^ bund similar lim- 



itations to the extent of the R-T fingers when radiative 
cooling was taken into account. 

The magnetic field is carried along by the plasma. Since 
the velocity field initially only has a radial component, the 
induction equation (9tB = V x (v x B)) implies that the 
toroidal field is swept up, while the radial field is just at- 
tenuated proportional to the square of the radius. In our 
models, the ejecta are initially set up using the same con- 
figuration for the magnetic field as is used for the CSM. 
As a result, the radial-field component is the dominant 
component in the ejecta. Starting at the R-T unstable re- 
gion at the contact discontinuity, a sizable toroidal field 
is present. The R-T instability induces a vg component, 
which in turn modifies the radial component of the mag- 
netic field. Along the R-T fingers, for Models B and C, 
the field is mostly radial with toroidal field components at 
the tips ar id bases of the fin gers, in agreement with results 
from e.g. lJun et al.l (|l995h . In Model A however, we do 
not see a significant radial component in the R-T unstable 
region, which may have to do with the small scale of the 
initial radial field. This case requires further investigation, 
allowing for a gradient in the third dimension. 

Models D, E, and F correspond to models A through 
C, but with an adiabatic index equal to 7 — 1.1. The com- 
pression ratio is analytically given by s = (7 + l)/(7 — 1) 
and therefore higher for a lower value of 7. As a conse- 
quence, the forward shock propagates more slowly into the 
CSM. At a time of 634 yr, the blast wave has only reached 
a distance of about 1.1 x 10"'^^ cm. For clarity purposes we 
show the SNR at a later stage, at an age of approximately 
1200 years. As a result of the higher compression ratio 
the forward and reverse shock are closer together and the 
R-T fingers almost reach the forward shock. The mag- 
netic field in models E and F again is obviously stretched 
around the R-T fingers, resulting in a dominantly radial 



field over the entire remnant. The magnetic field in model 
D, where the field was initially mainly toroidal, does not 
show a preferred radial component, even around the R-T 
fingers. 
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Figure 4: The SNR evolution is shown at a time of 1200 yrs after 
explosion. The initial magnetic field was mainly toroidal, and the 
adiabatic index 7 = 1.1 (Model D). The soft equation of state causes 
the relative proximity of forward to reverse shock, and the Rayleigh- 
Taylor fingers stretch almost out to the blast wave, seen best in the 
upper left quadrant, which shows the logarithm of the density. The 
lower left panel shows the logarithm of the absolute value of the 
radial magnetic field. The lower right part shows the logarithm of 
I I , and the upper right quadrant shows the logarithm of the ratio 
of the radial field component relative to the total field ( | Br | / 1 B | ) . 



4. Particle acceleration 

Observational evidence suggests that cosmic rays are 
efficiently accelerated at the forward shock, and p ossi- 
bly also at the reverse shock ( Helder and Vinkl . 2008f l. up 
to energies of 10^^ eV. The mechanism that is believed 
to be responsible for this is diffusive shock acceleration 
(DSA), also know as first-order Fermi acceleration. In 
this mechanism, the particle momentum increases every 
time the particle goes from upstream to downstream and 
back due to the randomization of the momentum vector 
by scattering in the rest frames of the up- and downstream 
flow. The simplest case of advection and diffusion of cos- 
mic rays in ID plane parallel geometry can be described 
by the following set of two sto chastic difference equations 
fiAchterberg and Krulls. . .l992. : ,van der Swaluw and Achterbere 
2004h : 



dV , 

du = —dt 

3dx 



(6) 
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Figure 5: Simulation of the SNR into CSM with a turbulent initial 
field and adiabatic index of 7 = 1.1 (Model E). The four panels 
represent the same as those in Figure |4] except for the lower right 
one, which represents the logarithm of \Bg\. The magnetic field is 
stretched radially around the R-T fingers, and reaches almost out to 
the forward shock. 



dx ^ Vdt + V2fidtC (7) 

The first equation represents the boost in energy (with 
u — la(j)/mc)), and the second one describes the random 
walk of the particle, with an advective term and a dif- 
fusive term. The diffusive term contains the variable ^, 
with (^i) = and = Sij, to represent the stochastic 

behaviour of the diffusion, of which the sum of many in- 
stances returns the distribution function of the relativistic 
particles. The diffusion of the particle is caused by inelas- 
tic scattering by a turbulent magnetic field, described by 
the diffusion parameter n. It depends on the energy of the 
particle, and the magnetic field strength, where diffusion 
parallel to the field may differ from the diffusion perpen- 
dicular to the field. In future work we will include this 
dependence, but for now we fix the diffusion parameter to 
a constant, in order to compare the spectrum to the an- 
alytical solution. The diffusive lengthscale XdiB = V'^Kdt 
needs to be larger than the shock thickness Xshock in order 
for the particles to feel the effects of DSA. The shock is 
spread over approximately 5 grid cells, which sets a lower 
limit to the diffusion parameter, depending on the resolu- 
tion. Additionally the advection length scale Xadv = Vdt 
should be smaller than the shock thickness, for the par- 
ticles to 'see' the shock. The condition that needs to be 
satisfied is 



AXadv < shock < ^^diff- 



(8) 



< 
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Figure 6: Simulation of SNR into CSM with an initial magnetic field 
parallel to the symmetry axis, i.e. vertically aligned (Model F). The 
four panels represent the same variables as in Figure[5] Similar to the 
case in model E, the magnetic field is radially dominant throughout 
the entire remnant, except for a small zone at the shock front itself. 



If these conditions are met, the resulting energy spec- 
trum should approach the analytical power-law solution 
of a strong shock, where slope q of the power-law is deter- 
mined by the shock co mpression ratio s = (7-I- l)/(7— 1). 
One has (jAchterberg and KruUsl . 1992 ): 



Fiu) oc p—r- oc p 
dp 



9+1 



with 



g=(,s + 2)/(s-l). 



(9) 



(10) 



In figure [7] we show the spectrum for the case where 
7 = 5/3, implicating s = 4 and q = 2. The coordinate sys- 
tem is set to slab geometry, such that no adiabatic losses 
occur. The mass and energy of the supernova explosion 
are respectively 2.5 Mq and 2.0 x 10^^ erg and the circum- 
stellar medium is taken to be homogeneous with a density 
of 1.64 X 10^^^ g cm^"^. Particles are introduced at the 
forward shock position, with a fixed initial energy, and the 
total number of particles increases linearly with time. 

The resulting spectrum for different times in the SNR 
evolution is plotted in Figure[71 For the higher energy par- 
ticles it takes longer to reach the equilibrium value of the 
spectrum that can be compared with the analytical solu- 
tion, and we can see the evolution of the spectrum towards 
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higher energies for later times. In the plane parallel geom- 
etry that we assume, there is no theoretical limit to the 
maximum energy the particles can attain. In spherical ge- 
ometry, in addition to the time during which particles have 
been accelerated, the maximum energy is limited by adi- 
abatic losses and t he size of the remnaiit (see e.g. iDrurv , 
198.4 lAchterberd . 120001 : iBell and Luce^ 1200 ih . DSA in 
spherical geometry will be treated and discussed in a fu- 
ture work. 
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Figure 7: Cosmic ray spectrum for 7 = 5/3, at different points in 
time of the SNR evolution. 



5. Conclusions and Discussion 

For an ideal mono-atomic gas that is not significantly 
influenced by the presence of cosmic rays, our models show 
that the radial magnetic field does not dominate in the 
outer region of the remnant, contrary to observations. For 
a softer equation of state, such as may be induced by cos- 
mic rays if particle acceleration is effective in the SNR, 
the field can actually become radial for the entire rem- 
nant, apart from a thin layer at the forward shock. This 
is seen in models E and F with the parallel and turbulent 
initial magnetic field topologies. Since cosmic rays are be- 
lieved to randomise the magnetic field near the shock in 
supernova remnants through streaming instabilities, the 
scenario where a magnetic field is initially turbulent and 
the equation of state is softer than 7 = 5/3 is a realis- 
tic possiblity ( Blondin and EllisonLl200l( ). Whether or not 
the equation of state can become even softer than that of a 
relativistic gas, where 7 = 4/3, depends on the slope of the 
spectrum. For a hard spectrum (g < 2), high energy par- 
ticles that escape from the system can carry away a sub- 
stantial fraction of the energy, thus lowering the effective 
adiabatic index. Observational evidence so far, favours a 
softer spectrum (q > 2), in which case the compression 
ratio would be hmited to up to about 7 for 7 — 4/3. Mod- 
els to explain the observed H.E.S.S. spectrum in 7-rays 



of the SNR RX J1713. 7-3946 bv 'Berez hko and Volkl (e.g. 
200i I2008I) that include nonlinear feedback between the 
CRs and the hydrodynamics, require a compression ratio 
of s = 6.3 in order to explain the observed spectrum. In 
such a case, the R-T fingers would not reach the forward 
shock. Although the toroidal field as used in models A 
and D is the favoured one for the magnetic field in stellar 
winds, this does not reproduce the radial field orientation 
from observations. It is possible however, that even in 
such a geometry, the cosmic ray streaming creates enough 
turbulence to effectively create a scenario like we use in 
model E. Even though further upstream the field may be 
toroidal, in the direct vicinity of the shock front and down- 
stream the field may be turbulent because of the cosmic 
rays, and consequently become mainly radial around the 
R-T fingers. The development of the radial field may be 
suppressed by the limitations of the 2.5D model, for which, 
because of the lack of a gradient in the 0-direction, the in- 
duction equation lacks a number of terms. This effect is 
present for all three models described here. Full 3D sim- 
ulations are needed to see if the results presented here 
are typical, including the possible difference in the devel- 
opment of R-T instabilities in 3D compared to the 2.5D 
results presented here. 

An alternative scenario to explain the radial magnetic_ 



field o rientation near the blast wave was proposed bv Zirakashvili and P 
( 2008h . They showed that the non-resonant streaming in- 
stability can amplify and radially stretch the magnetic field 
downstream of the shock front, when small density pertur- 
bations naturally create enough turbulence to trigger this 
effect in SNRs. The radial field can thus be a factor of 1.4 
larger than the component parallel to the shock front. 

Quantitatively, it is very complicated to make realistic 
predictions about values for the adiabatic index and level 
of turbulence. Since magnetic fields are addi tionally ampli- 
fied at the shock (e.g. lSell and Lucekl . l200ll ). the magnetic 
pressure may become an additional non-negligible compo- 
nent, again changing the compression ratio at the shock 
as well as the pressure profile. This was not taken into 
account in our models. In the test particle approach, the 
distribution of accelerated particles should follow a power 
law with a slope that is dependent on the adiabatic index. 
Preliminary results of our implementation of test particles 
and dynamical evolution with the MHD itself are promis- 
ing. Nonlinear feedback of the particle acceleration to the 
magnetic field and the equation of state are expected to 
be important. In future studies we hope to further evalu- 
ate the interaction and the resulting particle spectrum and 
observational features that it will produce. 
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